
#### APPENDIX C4: excluding outliers


# load, combine data
library(rio)
library(tidyverse)

# estimate models
library(brglm2)

# cluster standard errors
library(lmtest)
library(sandwich)

# generate tables
library(stargazer)



# LOAD DATA
glm_df <- import("glm_data.csv") %>%
  filter(iso3c!="TTO" & iso3c!="GUY" & iso3c!="SUR") %>%
  filter(!is.na(approval) & !is.na(opp1vote) & !is.na(number_protests))


## TABLE C6: ALL DOCS
brglm1 <- glm(doc ~ approval + opp1vote + number_protests +
                as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm1_r <- brglm1 %>% 
  coeftest(vcovHC(brglm1, type = 'HC0', cluster =  ~ iso3c))



brglm2 <- glm(doc ~ previous_doc + duck + partyage + polity2 + left_exec + any_election + minister_is_technocrat_mainstream + 
                as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm2_r <- brglm2 %>%
  coeftest(vcovHC(brglm2, type = 'HC0', cluster =  ~ iso3c))



brglm3 <- glm(doc ~ resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm3_r <- brglm3 %>% 
  coeftest(vcovHC(brglm3, type = 'HC0', cluster =  ~ iso3c))



brglm4 <- glm(doc ~ approval + opp1vote + number_protests +
                previous_doc + duck + partyage + polity2 + left_exec + any_election + minister_is_technocrat_mainstream + 
                as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm4_r <- brglm4 %>% 
  coeftest(vcovHC(brglm4, type = 'HC0', cluster =  ~ iso3c))



brglm5 <- glm(doc ~ approval + opp1vote + number_protests + 
                resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm5_r <- brglm5 %>%
  coeftest(vcovHC(brglm5, type = 'HC0', cluster = c('iso3c', 'Summary')))



brglm6 <- glm(doc ~ approval + opp1vote + number_protests + previous_doc + duck + partyage + polity2 + left_exec + any_election +
                minister_is_technocrat_mainstream + 
                resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm6_r <- brglm6 %>%
  coeftest(vcovHC(brglm6, type = 'HC0', cluster = c('iso3c', 'Summary')))

# generate tables
stargazer(brglm1_r,brglm2_r,brglm3_r,brglm4_r,brglm5_r,brglm6_r, #type="text",
          omit=c("iso3c","quarter","year"), no.space = T,
          covariate.labels = c("Executive approval (perc)","Opposition vote share (perc)","Number of protests",
                               "Previous document = 1","Lame duck = 1", "Party age",
                               "Democracy","Left executive = 1","Election quarter = 1",
                               "Mainstream minister = 1",
                               "Resource rents (perc)","Field discovery = 1","Oil price (USD, log)",
                               "GDP per capita (1,000 USD)","GDP growth (perc)","IMF program = 1", "Crisis = 1"), 
          dep.var.labels = "Any document")

# retrieve n, aic, ll, etc
stargazer(brglm1,brglm2,brglm3,brglm4,brglm5,brglm6, #type="text",
          dep.var.labels = "Any document")



## TABLE C7: FUNDS
brglm7 <- glm(fund_doc ~ approval + opp1vote + number_protests +
                as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm7_r <- brglm7 %>%
  coeftest(vcovHC(brglm7, type = 'HC0', cluster =  ~ iso3c))



brglm8 <- glm(fund_doc ~ previous_fund_doc + duck + partyage + polity2 + left_exec + any_election + minister_is_technocrat_mainstream + 
               as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
             family = binomial(link="logit"), 
             method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm8_r <- brglm8 %>% 
  coeftest(vcovHC(brglm8, type = 'HC0', cluster =  ~ iso3c))



brglm9 <- glm(fund_doc ~ resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = subset(glm_df, number_protests < 7))

brglm9_r <- brglm9 %>% 
  coeftest(vcovHC(brglm9, type = 'HC0', cluster =  ~ iso3c))



brglm10 <- glm(fund_doc ~  approval + opp1vote + number_protests +
                 previous_fund_doc + duck + partyage + polity2 + left_exec + any_election + minister_is_technocrat_mainstream + 
                 as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
               family = binomial(link="logit"), 
               method = "brglmFit",
               data = subset(glm_df, number_protests < 7))

brglm10_r <- brglm10 %>% 
  coeftest(vcovHC(brglm10, type = 'HC0', cluster =  ~ iso3c))



brglm11 <- glm(fund_doc ~ approval + opp1vote + number_protests + 
                 resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                 as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
               family = binomial(link="logit"), 
               method = "brglmFit",
               data = subset(glm_df, number_protests < 7))

brglm11_r <- brglm11 %>%
  coeftest(vcovHC(brglm11, type = 'HC0', cluster = c('iso3c', 'Summary')))



brglm12 <- glm(fund_doc ~ approval + opp1vote + number_protests + previous_fund_doc + duck + partyage + polity2 + left_exec + 
               any_election + minister_is_technocrat_mainstream + 
               resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
               as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
               family = binomial(link="logit"), 
               method = "brglmFit", maxit = 200, # need to run more iterations than the default (100) for the algorithm to converge
               data = subset(glm_df, number_protests < 7))

brglm12_r <- brglm12 %>%
  coeftest(vcovHC(brglm12, type = 'HC0', cluster = c('iso3c', 'Summary')))



# generate tables
stargazer(brglm7_r,brglm8_r,brglm9_r,brglm10_r,brglm11_r,brglm12_r, #type="text",
          omit=c("iso3c","quarter","year"), no.space = T,
          covariate.labels = c("Executive approval (perc)","Opposition vote share (perc)","Number of protests",
                               "Previous fund document = 1","Lame duck = 1", "Party age",
                               "Democracy","Left executive = 1","Election quarter = 1",
                               "Mainstream minister = 1",
                               "Resource rents (perc)","Field discovery = 1","Oil price (USD, log)",
                               "GDP per capita (1,000 USD)","GDP growth (perc)","IMF program = 1", "Crisis = 1"), 
          dep.var.labels = "Fund document")

# retrieve n, aic, ll, etc
stargazer(brglm7,brglm8,brglm9,brglm10,brglm11,brglm12, #type="text",
          omit=c("iso3c","quarter","year"), no.space = T,
          dep.var.labels = "Fund document")



